#ifndef A_DUAL_NUMBER_H_
#define A_DUAL_NUMBER_H_


/*!
    \page th_dual Dual numbers
    
    The concept of dual numbers is a useful expansion of a field of real
    numbers. Mathematically they are defined as elements of an algebra
    \f$ \Real[\eps] / \langle \eps^2 \rangle \f$. That means you can
    express any dual number as a polynomial \f$ x + \eps x' \f$, where
    \f$ x \f$ and \f$ x' \f$ are real numbers and \f$ \eps \f$ is
    a variable of the polynomial.
    
    Intuitively \f$ \eps \f$ can be imaged as a first-order small quantity,
    or as a real-like value which is less then any positive number.
    But you should keep in mind some nuances which can be essential
    if you takes this naive point of view. See \ref th_dual_extended
    for more information about them.
    
    \section th_dual_alg Algebraic operations
    There are some algebraic operations with the dual numbers, namely
    addition, subtraction and multiplication. Since any dual number is
    just a polynomial, you can do all these operations formally assuming
    \f$ \eps^2 = 0 \f$ (because the dual numbers form a factor ring by
    the ideal generated by the polynomial \f$ \eps^2 \f$):
    \f{align*}{
        (x + \eps x') + (y + \eps y') &= (x + y) + \eps (x' + y') \\
        (x + \eps x') - (y + \eps y') &= (x - y) + \eps (x' - y') \\
        (x + \eps x') \cdot (y + \eps y') &= x y + \eps (x' y + x y')
    \f}
    
    From the other hand, the dual numbers can be described as
    two-dimensional vectors over the field of real numbers. The canonical
    basis consist of two dual numbers \f$ 1 \f$ and \f$ \eps \f$ which are
    known as real and imaginary unit correspondingly, and the coefficients
    before them - as real and imaginary part of the given dual number:
    \f{align*}{
        \real(x + \eps x') &= x \\
        \imag(x + \eps x') &= x'
    \f}
    
    Notice the multiplication by a dual number \f$ a + \eps 0 \f$ with
    zero imaginary part is fully equivalent to the scaling by a real
    number \f$ a \f$:
    \f[
        (a + \eps 0) (x + \eps x') = a x + \eps (a x')
    \f]
    Therefore you can treat a dual number with zero imaginary part
    as the usual real number which is equal to its real part.
    
    \section th_dual_div Division
    For dual numbers division is defined as the inverse of multiplication.
    More accurate, the quotient of dual numbers \f$ a + \eps a' \f$ and
    \f$ b + \eps b' \f$ is a dual number \f$ x + \eps x' \f$ which
    satisfies the equation
    \f[
        a + \eps a' = (x + \eps x') (b + \eps b')
    \f]
    The equation can be transformed into the system of two linear
    equations:
    \f{align*}{
        a  &= x b \\
        a' &= x' b + x b'
    \f}
    
    In different cases that system has zero, one or infinity number
    of solutions:
    \li When \f$ b \neq 0 \f$ there is a single solution given by
    the formula:
    \f{align*}{
        x  &= a / b \\
        x' &= (a' / b) - (a / b) (b' / b)
    \f}
    \li When \f$ b = 0 \f$ and \f$ a \neq 0 \f$, the equations have
    no solution, and the result of division must contain infinity
    in the real part and infinity or NaN in the imaginary part.
    \li When \f$ a = b = 0 \f$, the system simplifies to a single
    equation:
    \f[
        a' = x b'
    \f]
    In this situation the solution is not unique, because one can
    substitute any number instead of \f$ x' \f$. So algebraically
    the quotient must contain NaNs in both real and imaginary parts.
    But in the extended mode the result is well-defined: it is a real
    number \f$ a' / b' \f$ (see corresponding \ref th_dual_extended
    "section" for discussion of the causes of this phenomenon).
    
    \section th_dual_function Functions of a dual argument
    Let \f$ F \f$ be analytic function of a single real argument. Then
    its value at the dual point \f$ x + \eps x' \f$ can be found via
    the Teylor series:
    \f[
        F(x + \eps x') =
        \sum_{n = 0}^{\infty} \frac{F^{(n)}(x)}{n!} (\eps x')^n =
        F(x) + \eps x' F'(x)
    \f]
    All terms with \f$ n \geq 2 \f$ are canceled because \f$ \eps \f$
    in a power of \f$ n \geq 2 \f$ is exactly zero.
    
    This idea can be generalized onto the wider class of functions.
    First of all, notice that in real numbers evaluation of a function
    at the specified point is a linear functional on the space of
    functions. Analogically any dual number \f$ x + \eps x' \f$ generates
    a pair of linear functionals (which will be noticed as \f$ \pi \f$
    and \f$ \delta \f$) for evaluation a real and imaginary part of
    the result:
    \f[
        F(x + \eps x') =
        \langle \pi | F \rangle +
        \eps \langle \delta | F \rangle
    \f]
    The theory of functions of dual argument is based on the three axioms:
    \li Functional \f$ \pi \f$ is multiplicative:
    \f$
        \langle \pi | FG \rangle =
        \langle \pi | F \rangle
        \langle \pi | G \rangle
    \f$
    \li Functional \f$ \delta \f$ submits to the Leibnitz law:
    \f$
        \langle \delta | FG \rangle =
        \langle \delta | F \rangle \langle \pi | G \rangle +
        \langle \pi | F \rangle \langle \delta | G \rangle
    \f$
    \li 
    
    
    Analogically, if \f$ F \f$ is a multivariate function, then
    the imaginary part of its value at \f$ x + \eps x' \f$ (where
    \f$ x \f$ and \f$ x' \f$ are vectors) contains directional derivative
    along the vector \f$ x' \f$.
    
    In the file dual_functions.h there are implementation of common-used
    mathematical functions for dual numbers. It's a good idea to combine
    them into expressions for automatical calculation of the derivative.
    
    \section th_dual_singleside Single-side mode
    There are many functions which are not differentiable in some points,
    but have a single-side derivative there. For example, the function of
    absolute value \f$ z = |x| \f$ at zero has no derivative, but it has
    both right-side and left-side derivatives, whose are equal to
    +1 and -1 correspondingly. In the common case, single-side derivatives
    of a function \f$ F \f$ are defined by the following:
    \f{align*}{
        F^{+}(x) &= \lim_{\nu \to 0+} \frac{F(x + \nu) - F(x)}{\nu} \\
        F^{-}(x) &= \lim_{\nu \to 0-} \frac{F(x + \nu) - F(x)}{\nu}
    \f}
    
    For the function of several arguments the notions of the left-side and
    right-side derivatives are meaningless (because in the multidimensional
    space there is no "left" and "right"), but the idea of single-side
    derivation still works good. So, we can tell about the single-side
    directional derivative along the vector \f$ x' \f$:
    \f[
        F^{+}(x) = \lim_{\nu \to 0+} \frac{F(x + \nu x') - F(x)}{\nu}
    \f]
    
    An important example of a single-side derivation of a multivariate
    function is a derivative of the minimum function \f$ \min(x, y) \f$
    at the point \f$ x = y \f$:
    \f{align*}{
        F^{+}(x, x) &=
        \lim_{\nu \to 0+}
            \frac{\min(x + \nu x', x + \nu y') - \min(x, x)}{\nu} =\\&=
        \lim_{\nu \to 0+}
            \frac{x + \nu \min(x', y') - x}{\nu} =\\&=
        \min(x', y')
    \f}
    The minimum function plays a major role on the ordering of dual
    numbers. Dual numbers equpped with the minimum function form
    a semilattice which induces a total order. Namely, the dual
    numbers \f$ \alpha \f$ and \f$ \beta \f$ are ordered as
    \f$ \alpha \leq \beta \f$ if and only if \f$ \min(\alpha, \beta) =
    \alpha \f$. See dual_number::compare() for more information.
    
    Implementation of dual numbers in the Alib supposes to make a deal
    with both single-side and double-side derivatives. The flag
    dual_traits::single_side_mode controls the mode of derivation:
    if it is \c true, imaginary part will store a single-side derivative,
    otherwise it will store a double-side derivative or NaN if
    the function has a fracture at the specified point.
    
    \section th_dual_extended Extended mode
    In this situation, there are two competing approaches. The first one
    is based 
*/


//! A class that represents dual number in Cartesian form.
/*!
    Dual number (or hypercomplex number of parabolic type) is an ordered
    pair of real numbers. The first one is called real part, and the
    second one is imaginary part. Any dual number can be expressed
    in form \f$ a + \eps b \f$, where \f$ a \f$ and \f$ b \f$ are
    real and imaginary part of the dual number, and \f$ \eps \f$
    is a formal symbol, which is known as imaginary unit.
    
    All operations with dual numbers can be built formally (algebraically)
    as operations with a polynomials \f$ a + \eps b \f$ with real
    coefficients, where
    \f[
        \boxed{\eps^2 = 0}
    \f]
    This formula is a main equality of dual numbers. That's why it is
    marked with a bounded rectangle. Please pay special attention on it.
    
    On your practice you can look to the imaginary unit like to a small
    value of the first order. In other words, \f$ \eps \f$ is a real-like
    quantity, which is less then any real positive number, and its square
    is indistinguishable of zero.
    
    \tparam scalar A type of real numbers, which must supply a full
    complement of algebraic and ordering operators.
    
    \tparam traits Collection of special properties of real and dual
    numbers. See \ref dual_traits class for more information.
    
    \todo Rewrite this class as derived of a common class of hypercomplex
    numbers.
*/
template
<
    class scalar,
    class traits = dual_traits<scalar>
>
class dual_number
{
private:
    scalar re, im;

public:
    //! Default constructor.
    /*!
        Construct dual number with no initialization.
    */
    dual_number()
    {
    }
    
    
    //! General constructor.
    /*!
        Construct dual number \f$ x + \eps y \f$ from specified real
        and imaginary parts \a x and \a y. You can omit imaginary part,
        because it is zero by default. In other words, this constructor
        builds real number if you don't specify an imaginary part.
        
        \note The constructor is not declared as \c explicit, so
        implicit conversion from real to dual number is possible,
        but not recommended.
    */
    dual_number( scalar const & x, scalar const & y = 0 )
    {
        this->re = x;
        this->im = y;
    }
    
    
    //! Extract the real part of \c this dual number.
    /*!
        \math Mathematically this operation is given by the formula:
        \f[
            \real(a + \eps b) = a
        \f]
        
        \sa imag().
    */
    scalar const & real() const
    {
        return re;
    }
    
    
    //! Extract the imaginary part of \c this dual number.
    /*!
        \math Mathematically this operation is given by the formula:
        \f[
            \imag(a + \eps b) = b
        \f]
        
        \sa real().
    */
    scalar const & imag() const
    {
        return im;
    }
    
    
    //! Negate \c this dual number and store result to itself.
    /*!
        \math Mathematically this operation is given by the formula:
        \f[
            -(a + \eps b) = (-a) + \eps (-b)
        \f]
        
        \return Reference to \c this instance, storing a result.
        
        \sa operator-().
    */
    dual_number & negate()
    {
        re = -re;
        im = -im;
        return *this;
    }
    
    
    //! Conjugate \c this dual number and store result to itself.
    /*!
        \math Mathematically this operation is given by the formula:
        \f[
            (a + \eps b)^{\dag} = a - \eps b
        \f]
        
        \return Reference to \c this instance, storing a result.
        
        \sa operator~().
    */
    dual_number & conj()
    {
        im = -im;
        return *this;
    }
    
    
    //! Multiplication of \c this dual number by a real number.
    /*! 
        The function computes product of \c this dual number
        and a specified real number \a x using the formula below,
        places result in \c this instance, and returns reference to it.
        
        \math Dual numbers form a linear space over the field of real
        numbers. Therefore multiplication of dual number
        \f$ a + \eps b \f$ and a real number \f$ x \f$ is a per-component
        operation:
        \f[
            (a + \eps b) \cdot x = ax + \eps (bx)
        \f]
        
        \sa
        \li \mul_dual_dual      "mul(dual_number, dual_number)",
        \li \op_mul_scalar      "operator*=(scalar)",
        \li \op_mul_dual_scalar "operator*(dual_number, scalar)",
        \li \op_mul_scalar_dual "operator*(scalar, dual_number)".
    */
    dual_number & mul( scalar const & x )
    {
        re *= x;
        im *= x;
        return *this;
    }
    
    
    //! Division of \c this dual number by a real number.
    /*!
        The function computes fraction of \c this dual number
        and a specified real number \a x using the formula below,
        places result in \c this instance, and returns reference to it.
        
        \math Dual numbers form a linear space over the field of real
        numbers. Division of dual number \f$ a + \eps b \f$ by a real
        number \f$ x \f$ fully equivalent to \mul_scalar "multiplication"
        by \f$ 1 / x \f$:
        \f[
            \frac{a + \eps b}{x} = \frac{a}{x} + \eps \frac{b}{x}
        \f]
        
        \sa
        \li \div_dual_dual      "div(dual_number, dual_number)",
        \li \op_div_scalar      "operator/=(scalar)",
        \li \op_div_dual_scalar "operator/(dual_number, scalar)",
    */
    dual_number & div( scalar const & x )
    {
        re /= x;
        im /= x;
        return *this;
    }
    
    
    //! Addition of two dual numbers, storing result in \c this.
    /*!
        The function computes sum of specified dual numbers \a x and \a y
        using the formula below, places result in \c this instance,
        and returns reference to it.
        
        \math The sum of two dual numbers is defined algebraically
        via terms regrouping:
        \f[
            (a + \eps b) + (x + \eps y) = (a + x) + \eps (b + y)
        \f]
        
        \sa
        \li \op_add_dual_dual "operator+(dual_number, dual_number)",
        \li \op_add_dual      "operator+(dual_number)".
    */
    dual_number & add( dual_number const & x, dual_number const & y )
    {
        re = x.re + y.re;
        im = x.im + y.im;
        return *this;
    }


    //! Subtraction of two dual numbers, storing result in \c this.
    /*!
        The function computes difference of specified dual numbers \a x
        and \a y using the formula below, places result in \c this
        instance, and returns reference to it.
        
        \math The difference of two dual numbers is defined pure
        algebraically via terms regrouping:
        \f[
            (a + \eps b) - (x + \eps y) = (a - x) + \eps (b - y)
        \f]
        
        \sa
        \li \op_sub_dual_dual "operator-(dual_number, dual_number)",
        \li \op_sub_dual      "operator-=(dual_number)".
    */
    dual_number & sub( dual_number const & x, dual_number const & y )
    {
        re = x.re - y.re;
        im = x.im - y.im;
        return *this;
    }
    
    
    //! Multiplication of two dual numbers, storing result in \c this.
    /*!
        The function computes product of specified dual numbers \a x
        and \a y using the formula below, places result in \c this
        instance, and returns reference to it.
        
        \math The product of two dual numbers is defined algebraically
        considering an equality \f$ \eps^2 = 0 \f$:
        \f[
            (a + \eps b) - (x + \eps y) = (a - x) + \eps (b - y)
        \f]
        Note that this definition has no contradiction with \mul_scalar
        "multiplication by real number". It means, after substitution
        \f$ y = 0 \f$ to the last formula you get the formula for scalar
        multiplication, and the following lines causes the same result:
        \code
            dual.mul(0.3);                      // dual *= 0.3
            dual.mul(dual, dual_number(0.3));   // dual *= (0.3 + 0 * eps)
        \endcode
        
        \note All essential code of this method is concentrated in
        \c traits::multiply(). You can reimplement that function
        if you need fine tuning of behavior. For instance, overriding
        \c traits::multiply() with SSE features allows to improve
        performance or accuracy (see \ref dual_traits::multiply()
        as an example).
        
        \sa
        \li \op_mul_dual_dual "operator*(dual_number, dual_number)",
        \li \op_mul_dual      "operator*=(dual_number)".
    */
    dual_number & mul( dual_number const & x, dual_number const & y )
    {
        traits::multiply(&this, x, y);
        return *this;
    }


    //! Division of two dual numbers, storing result in \c this.
    /*!
        The function computes fraction of specified dual numbers \a x
        and \a y using the formula below, places result in \c this
        instance, and returns reference to it.
        
        \math Let's first describe fraction of real unit and an arbitrary
        dual number, which is known as inversion. The inverse of the dual
        number \f$ a + \eps b \f$ is a dual number \f$ x + \eps y \f$,
        whose product with \f$ a + \eps b \f$ is a simply one:
        \f[
            (a + \eps b)(x + \eps y) =
            (ax) + \eps (bx + ay) = 1
        \f]
        From the last equality, we get:
        \f[
            x = \frac{1}{a} \qquad
            y = -\frac{b}{a^2}
        \f]
        Note that the inverse of \f$ a + \eps b \f$ can be written as
        \f[
            (a + \eps b)^{-1} = \frac{(a + \eps b)^{\dag}}{a^2}
        \f]
        The formula shows there are non-zero dual numbers, which
        cannot be inverted. All of them has zero real part.
        
        \cont Now we define fraction of two dual numbers as a product of
        dividend and inversion of divisor:
        \f[
            \frac{a + \eps b}{x + \eps y} =
            \frac{(a + \eps b) (x - \eps y)}{x^2} =
            \frac{a}{x} + \eps \frac{b - a y / x}{x}
        \f]
        Note that this definition has no contradiction with \div_scalar
        "division by real number". It means, after substitution
        \f$ y = 0 \f$ to the last formula you get the correct formula
        for scalar division, and the following lines causes the same
        result:
        \code
            dual.div(0.3);                      // dual /= 0.3
            dual.div(dual, dual_number(0.3));   // dual /= (0.3 + 0 * eps)
        \endcode
        
        \note All essential code of this method is concentrated in
        \c traits::divide(). You can reimplement that function
        if you need fine tuning of behavior. For instance, overriding
        \c traits::divide() with SSE features allows to improve
        performance or accuracy (see \ref dual_traits::divide()
        as an example).
        
        \sa
        \li \op_div_dual_dual "operator/(dual_number, dual_number)",
        \li \op_div_dual      "operator/=(dual_number)".
    */
    dual_number & div( dual_number const & x, dual_number const & y )
    {
        traits::divide(&this, x, y);
        return *this;
    }
    
    
    //! Compare specified dual numbers.
    /*!
        \math Unlike complex and split numbers, dual numbers can
        be well-ordered. This ordering is connected with representation
        an imaginary unit as a small quantity of the first order. Dual
        number \f$ a + \eps b \f$ is less then \f$ x + \eps y \f$,
        if \f$ a < x \f$, or if its real parts are equal and \f$ b < y \f$.
        As you see, this is regular lexicographical ordering of pairs
        of real numbers.
    */
    static typename traits::compare_result
        compare( dual_number const & x, dual_number const & y )
    {
        //! \todo Implement and write documentation.
    }
    
    
    //! Unary operator +, does nothing.
    /*!
        The operator simply returns a copy of \c this instance.
    */
    dual_number operator+() const
    {
        return dual_number(*this);
    }
    
    
    //! Unary operator - (negation).
    /*!
        The operator returns a dual number, which is negative of \c this
        and keeps \c this without changing.
        
        \sa negate().
    */
    dual_number operator-() const
    {
        return dual_number(*this).negate();
    }
    
    
    //! Unary operator ~ (conjunction).
    /*!
        The operator returns a dual number, which is negative of \c this
        and keeps \c this without changing.
        
        \sa conj().
    */
    dual_number operator~() const
    {
        return dual_number(*this).conj();
    }
    
    dual_number const & operator=( scalar const & x )
    {
        re = x;
        im = 0;
    }
    
    
    dual_number const & operator+=( scalar const & x )
    {
        return add(*this, dual_number(x));
    }
    
    dual_number const & operator-=( scalar const & x )
    {
        return sub(*this, dual_number(x));
    }
    
    dual_number const & operator*=( scalar const & x )
    {
        return mul(x);
    }
    
    dual_number const & operator/=( scalar const & x )
    {
        return div(x);
    }


    //! Addition operator
    dual_number const & operator+=( dual_number const & x )
    {
        return add(*this, x);
    }
    
    dual_number const & operator-=( dual_number const & x )
    {
        return sub(*this, x);
    }
    
    dual_number const & operator*=( dual_number const & x )
    {
        return mul(*this, x);
    }
    
    dual_number const & operator/=( dual_number const & x )
    {
        return div(*this, x);
    }
};


template< class scalar, class traits >
dual_number<scalar, traits> const &
    operator+( dual_number<scalar, traits> const & x, scalar const & y )
{
    return dual_number<scalar, traits>()
        .add(x, dual_number<scalar, traits>(y, 0));
}


#endif /* A_DUAL_NUMBER_H_ */
